function [p,w,E_i,E_r,LMat, VMat]= PriceWages(pm,ze,zr)

seedha = zr>ze;

thresh = 100;
p =ones(pm.I,1);
P = exp(sum(pm.kappa.*log(p./pm.kappa)));
p = p./P;


Ee_i= zeros(pm.I,1);
for r = 1:pm.I
    Ee_i(r)= quadgk(@(e) e.^(pm.theta(r)./(pm.theta(r)-pm.rho)).*lognpdf(e,0,pm.se),-Inf,Inf);
end

w=10;
tol = 10; iter = 1; conv = 0; TolMat = zeros(150,1);
%pm.bp = pm.ratio.^(1./pm.nu).*pm.bc;


while tol>pm.mintolp && iter<200 %|| (conv>=pm.I-1 && iter>100)
%        wf_zr = wf(pm,zr,w);
        %wages and implied fixed costs 
        E_i = zeros(pm.I,1);
        E_r = zeros(pm.I,1);
        for r = 1:pm.I
            E_i(r) = integral(@(e) pi_i(pm,ze(r),e,p(r),w,r).*lognpdf(e,0,pm.se),-Inf,Inf)./p(r);
            E_r(r) = integral(@(e) (pi_f(pm,zr(r),e,p(r),w,pm.t(r),r) - pi_i(pm,zr(r),e,p(r),w,r)) ...
                    .*lognpdf(e,0,pm.se),-Inf,Inf)./p(r);
        end

        %Total profits in the formal and informal sectors 
        L_i = zeros(pm.I,1);
        L_f = zeros(pm.I,1);
        L_c = zeros(pm.I,1);
        L_p = zeros(pm.I,1);
        wfL_f = zeros(pm.I,1);
        Y_i = zeros(pm.I,1);
        Y_f = zeros(pm.I,1);
        Pi_i = zeros(pm.I,1);
        Pi_f = zeros(pm.I,1);
        
      
        for r = 1:pm.I
            L_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  l_i(pm,x,e,p(r),w,r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            L_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  l_f(pm,x,e,p(r),w,pm.t(r),r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            wfL_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x)  ...
                    wf(pm,x,e,w,r).*l_f(pm,x,e,p(r),w,pm.t(r),r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            Y_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) x.*e.* l_i(pm,x,e,p(r),w,r).^pm.rho...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            Y_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) x.*e.*  l_f(pm,x,e,p(r),w,pm.t(r),r).^pm.rho...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            Pi_i(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) (pi_i(pm,x,e,p(r),w,r) - p(r).*E_i(r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf);
            
            Pi_f(r) = pm.M(r).*seedha(r).*integral(@(e) integral(@(x) (pi_f(pm,x,e,p(r),w,pm.t(r),r) - p(r).*(E_i(r)+E_r(r)))...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);

            L_c(r) = integral(@(e) integral(@(x) ...
                    (w.*pm.bc(r)./((x.*e).^pm.phi.*wf(pm,x,e,w,r))).^(-pm.nu) .*L_f(r) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),Inf,'ArrayValued',true),-Inf,Inf);

            L_p(r) = integral(@(e) integral(@(x) ...
                    (w.*pm.bp(r)./wf(pm,x,e,w,r)).^(-pm.nu) .*L_f(r) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf);
        end


        
        %Fixed costs and Taxes paid
        FC = p.*E_i.*pm.M.*(1-logncdf(min(ze,zr),0,pm.sigma)) + p.*E_r.*pm.M.*(1-logncdf(zr,0,pm.sigma));
        Tpaid = pm.t.*p.*Y_f;

%        income = w.*sum(L_i + pm.bc.*L_c + pm.bp.*L_p,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
        income = sum(w.*L_i + wfL_f,'omitnan') + sum(Pi_i + Pi_f + Tpaid,'omitnan');
        pnew = (pm.kappa.*income + FC)./ (Y_i+Y_f);
        pnew(pnew>10)=10;
        Pnew = exp(sum(pm.kappa.*log(pnew./pm.kappa)));
        pnew = pnew./Pnew;

        frac = sum(L_i+ L_c + L_p,'omitnan')./pm.L;
        tolp = abs((p-pnew))./p;
        tolw = abs(frac-1);
        a= 0.7;
        p = a.*p +(1-a).*pnew;
        w = a.*w +(1-a).*frac.*w;
        
%        p(1) = pm.alpha(1).*exp(-1./pm.alpha(1) .*sum(pm.alpha(2:pm.I).*log(p(2:pm.I)./pm.alpha(2:pm.I))));
        tol = max(max(tolp),tolw);
        conv = sum(tolp<=pm.mintolp)+sum(tolw<=pm.mintolp);
        %fprintf('Iter = %d; Converged = %d out of %d; tol = %0.4f \n',iter,conv,pm.I+1,tol);
        
        %End the algorithm if the guess is not working
        TolMat(iter) = tol;

%        if iter>=1+thresh || sum(TolMat(iter-thresh+1:iter))==thresh
%           iter= 1e3-1;
%            tol = 0;
%        fprintf('Not converging \n');
%        end
%        [tolw;tolp]
        iter = iter+1;
end

LMat = [L_i L_c+L_p];
%fprintf(', Guess Completed in %d and tol = %0.4f \n',iter, tol);
%fprintf('Guess Completed in %d \n',iter);

varz_i = zeros(pm.I,1);
varz_f = zeros(pm.I,1);
for r = 1:pm.I
            lbar_i = integral(@(e) integral(@(x)  log(l_i(pm,x,e,p(r),w,r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma));

            lbar_f = integral(@(e) integral(@(x)  log(l_f(pm,x,e,p(r),w,pm.t(r),r)) ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),5.*zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(1 - logncdf(zr(r),0,pm.sigma));

            %variance of x
            varz_i(r) = integral(@(e) integral(@(x)  (log(l_i(pm,x,e,p(r),w,r)) - lbar_i).^2 ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma));

            varz_f(r) = integral(@(e) integral(@(x)  (log(l_f(pm,x,e,p(r),w,pm.t(r),r)) - lbar_f).^2 ...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),5.*zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(1 - logncdf(zr(r),0,pm.sigma));
end
VMat = [varz_i; varz_f]; 
